Novelty metrics

Estimation of climate velocity (magnitude and direction) for the Phytoclimate of multiple growth forms
Author

Alejandro Ordonez

Published

March 22, 2023

The setup.

Here, I will be exploring how the novelty metrics in (Ordonez, Williams, and Svenning 2016) can be used in (Conradi et al. 2020) work on Phyto-climates that builds on his work on operationalizing the definition of the biome for global change research.

Input information.

Here using will use the maps of growth form (GF) suitability for 14 GFs:

Acronym Name
TE Evergreen trees
TDdry Drought-deciduous trees
TDcold Cold-deciduous trees
TN Needleleaf trees
ShE Evergreen shrubs
ShDdry Drought-deciduous shrubs
ShDcold Cold-deciduous shrubs
H Herbs
Geo Geophytes
Thero Therophytes
GC3 C3 grasses
GC4 C4 grasses
Suc Succulents
Clim Climbers

Compositional Novelty as a metric of Ecological Novelty.

The novelty metrics in (Ordonez, Williams, and Svenning 2016) focus on measuring three different mechanisms by which ecological novelty might emerge:

This mechanism is based on the idea that as environmental changes happen the composition of taxa on a site will change change, and therefore you would like to assess of this change means that this assemblage is new when compared to past assemblages.

To measure this, I focus on the composition rearrangement of a location due to environmental changes. In those situations where current communities that are compositionally (significantly) different to those found in the past, it can be considered that these assemblages are “novel”. Conceptual this idea links to (Williams, Jackson, and Kutzbach 2007) and (Williams and Jackson 2007) discussions on Novel climates and no-analog communities.

Core to measuring novelty is the direction of the contrasts. So Measuring if current conditions are novel means contrasting Each current condition, to ALL past cells. Therefore, to start I load the suitability per-growth form under the current environmental conditions (here defined as 1950 values, Figure 1):

Code
## Make suitability raster for all growth forms for the present time step [1950] - 45th slot in the data list 
Maps0kaBPlist <- lapply(dimnames(pml[[45]])[[2]],
                       function(i){
                         rasterize(x = as.matrix(xy), # points as a matrix
                                   y = cr.ea, #  template raster
                                   values = pml[[45]][, i] # Values to aggregate
                                   )
                         })
# Turn the list into a multi-layer SpatRaster
Maps0kaBP <- do.call("c",Maps0kaBPlist)
# Add the GF names to each layer
names(Maps0kaBP) <- NamesDtFrm$Name[match(dimnames(pml[[1]])[[2]],NamesDtFrm$Acro2)]

# Plot the Suitability Raster
levelplot(raster::stack(Maps0kaBP), # level plot only works with raster files
          at=seq(0,0.6,length.out=100), # Set the zlim
          scales = list(draw=FALSE), # To remove the Latlong
          main=list("Suitability per growth form [1950]",side=1,line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(col=rev(hcl.colors(100, palette = "RdYlBu"))) # add a legend
          ) + 
# Add the country outlines
layer(sp.lines(wrld_simpl)
      )

## Generate a data.frame of cells with data for present conditions
Maps0kaBPTble <- values(Maps0kaBP,
                        dataframe = TRUE,
                        na.rm = TRUE)

Figure 1: Present (1950) suitability per-growth form.

After this, I load the suitability per-growth form under past environmental conditions (Figure 2), here defined as the Last Glacial Maximum (21kaBP) values:

Code
## Make suitability raster for all growth forms for the first time step [21kaBP?]
Maps21kaBPlist <- lapply(dimnames(pml[[1]])[[2]],
                       function(i){
                         rasterize(x = as.matrix(xy), # points as a matrix
                                   y = cr.ea, #  template raster
                                   values = pml[[1]][, i] # Values to aggregate
                                   )
                         })
# Turn the list into a multi-layer SpatRaster
Maps21kaBP <- do.call("c",Maps21kaBPlist)
# Add the GF names to each layer
names(Maps21kaBP) <- NamesDtFrm$Name[match(dimnames(pml[[1]])[[2]],NamesDtFrm$Acro2)]
# Plot the map
# Plot the Suitability Raster
levelplot(raster::stack(Maps21kaBP), # level plot only works with raster files
          at = seq(0,0.6,length.out=100), # Add limits
          scales = list(draw=FALSE), # To remove the Latlong
          main=list("Suitability per growth form [21kaBP]",side=1,line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(col=rev(hcl.colors(100, palette = "RdYlBu")))  # add a legend
          ) + 
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

## Generate a data.frame of cells with data for past conditions
Maps21kaBPTble<- values(Maps21kaBP,
                        dataframe = TRUE,
                        na.rm = TRUE)

Figure 2: Past (21kaBP) Suitability per-growth form

The values in the two figures above are defined as the proportion of the species within a growth form for which the environmental conditions at 50x50km a grid-cell are considered suitable at a given period (here 1950 and 21kaBP).

This Suitability is based on a Eco-physiological species distribution Model (INSERT NAME HERE).

Now with the adequate temporal data (i.e., Current & Past conditions), I estimate novelty. As stated above, the approach used to estimate novelty (EACH current conditions vs. ALL past cells approach used) means that for each current cell, I will have a large number of possible contrasts based on an adequate distance metric (e.g., (Squared) chord-distance or \(\chi^2\)-distance for composition data; Euclidean distance for uncorrelated environmental data, or Mahalanobis distance for correlated environmental data).

Given that the data I have for growth-form suitability is the proportion of species within a group for which the evaluated cell has suitable conditions (similar to proportion of species) I will use the Min Chord distance as suggested by (Simpson 2007), (Overpeck, Webb, and Prentice 1985), and (Gavin et al. 2003).

The chord distance between samples \(j\) and \(k\), \(d_{jk}\), is:

\[d_{jk} = \sqrt{\sum_{k=1}^{m}( x_{jk}^{0.5} - x_{ik}^{0.5})^2}\] Where \(x_{ij}\) is the proportion of growth from \(i\) in sample \(k\). Note that other dissimilarity metrics could be used, and they are implemented in the anaolgue package.

With all the pairwise distances estimated, a technique used to estimate site novelty (as in (Williams, Jackson, and Kutzbach 2007), (Williams and Jackson 2007), (Ordonez et al. 2014), and (Ordonez, Williams, and Svenning 2016)) is to retain the minimum dissimilarity value of the contrast between the target assemblage (here 1950’s sites) and all sites in the reference period (here 21kaBP). This technique is similar to the analogue matching approach in paleoecology ((Overpeck, Webb, and Prentice 1985) and (Flower, Juggins, and Battarbee 1997)).

Now using the Min Chord distance (as implemented in anaolgue), I estimate (using a parallel computing approach) the novelty of each cell in Maps0kaBP the present growth from suitability measured based on 1950 climates, 3 when compared to the Last Glacial Maximum growth from suitability as presented in Maps21kaBP:

Code
### Calculate the Min Chord distance for each current cell in Maps0kaBP to all past cells in Maps21kaBP
if(!"CordD_min.tif"%in%dir("./Data/Novelty/")){
  a<-Sys.time()
# Create a virtual cluster
  sfInit(cpus=10,parallel=TRUE)
## Export packages
  sfLibrary(analogue)
## Export Data
  sfExport("Maps0kaBPTble")
  sfExport("Maps21kaBPTble")

# Calculate the Squared Chord distance for each current cell in Maps0kaBP to all past cells in Maps21kaBP
  CordDistSumList <- sfLapply(1:dim(Maps0kaBPTble)[1],#(x<-1)
                              function(x, DistMet = "chord"){
                                DistCalc <- analogue::distance(x = Maps0kaBPTble[x,],
                                                               y = Maps21kaBPTble,
                                                               method = DistMet,
                                                               double.zero = TRUE)
                                
                                min(DistCalc,na.rm=T)
                              })
  sfStop();gc()
  CordDistSum <- data.frame(Dist = do.call("c",CordDistSumList),
                            CellID = as.numeric(rownames(Maps0kaBPTble)))
  write.csv(CordDistSum,"./Data/Novelty/CordDist_min.csv")
# Turn the Min Chord distance per cell into a raster
#CordDistSum <- read.csv("./Data/CordDist_min.csv") # Load data
  CordDistRast <- rast(Maps0kaBP[[1]]) # create the empty raster
  values(CordDistRast) <- NA #| fill with empty data
  names(CordDistRast) <- "minCordDist" # change layer name
  values(CordDistRast)[CordDistSum$CellID] <- CordDistSum$Dist # add the ChordD.min data
# Save the Raster file
  writeRaster(CordDistRast,
              "./Data/Novelty/CordD_min.tif",
              overwrite=TRUE)
  Sys.time()-a # Should clock about 3 to 5 minutes
}

# Load the CordDistRast raster if you don't run the Min Chord distance per cell estimation
if(!"CordDistRast"%in%ls()){
  CordDistRast <- rast("./Data/Novelty/CordD_min.tif")
}
# plot the minimum chord distance
plot(CordDistRast,
     main = "Min chord distance",
     col = rev(hcl.colors(100,"RdYlBu")))
plot(wrld_simpl, add=T)

Figure 3: Min Chord distance for each current cell in Maps0kaBP to all past cells in Maps21kaBP.

I also do this estimation using a Chi-Squared distance (Figure 4):

Code
### Calculate the Min Chi-Sqred distance for each current cell in Maps0kaBP to all past cells in Maps21kaBP
if(!"CordD_min.tif"%in%dir("./Data/Novelty/")){
  a<-Sys.time()
# Create a virtual cluster
  sfInit(cpus=10,parallel=TRUE)
# Export packages
  sfLibrary(analogue)
# Export Data
  sfExport("Maps0kaBPTble")
  sfExport("Maps21kaBPTble")
  
# Calculate the Squared Chord distance for each current cell in Maps0kaBP to all past cells in Maps21kaBP
  ChiDistSumList <- sfLapply(1:dim(Maps0kaBPTble)[1],#(x<-1)
                              function(x, DistMet = "chi.square"){
                                DistCalc <- analogue::distance(x = Maps0kaBPTble[x,],
                                                               y = Maps21kaBPTble,
                                                               method = DistMet,
                                                               double.zero = TRUE)
                                
                                min(DistCalc,na.rm=T)
                              })
  sfStop()
  ChiDistSum <- data.frame(Dist = do.call("c",ChiDistSumList),
                            CellID = as.numeric(rownames(Maps0kaBPTble)))
  write.csv(ChiDistSum,"./Data/Novelty/ChiDist_min.csv")
  
# Turn the Min Chord distance per cell into a raster
# ChiDistSum <- read.csv("./Data/ChiDist_min.csv") # Load data
  ChiDistRast <- rast(Maps0kaBP[[1]]) # create the empty raster
  values(ChiDistRast) <- NA #| fill with empty data
  names(ChiDistRast) <- "minChiDist" # change layer name
  values(ChiDistRast)[ChiDistSum$CellID] <- ChiDistSum$Dist # add the ChordD.min data
# Save the Raster file
  writeRaster(ChiDistRast,
              "./Data/Novelty/ChiD_min.tif",
              overwrite=TRUE)
  Sys.time()-a # Should clock about 3 to 5 minutes
}

# Load the CordDistRast raster if you don't run the Min Chord distance per cell estimation
if(!"ChiDistRast"%in%ls()){
  ChiDistRast <- rast("./Data/Novelty/ChiD_min.tif")
}
# plot the minimum chord distance
plot(ChiDistRast,
     main = "Min Chi-Sqr distance",
     col = rev(hcl.colors(100,"RdYlBu")))
plot(wrld_simpl, add=T)

Figure 4: Min Chi-Sqred distance for each current cell in Maps0kaBP to all past cells in Maps21kaBP.

The Next step in producing a novelty map is defining the a suitable dissimilarity cut-off to determine if the composition difference (i.e. the minimum distance) means a current assemblages ar novel in contrast to those in the past.

One way to do this, as suggested by (Simpson 2007) in the manual for anaolgue, is to use Monte Carlo simulation to determine a dissimilarity threshold that is unlikely to have occurred by chance. For this, two samples are drawn, at random, from the training set (i.e., the modern sample) and the dissimilarity between these two samples is recorded. This process is repeated many times to generate a randomization distribution of dissimilarity values expected by random comparison of samples. The dissimilarity value that achieves at a significance level of 0.01 can be determined by selecting the 0.01 probability percentile of the randomization distribution (the 1st percentile). Is important to note that to define this value at a 0.01 significance level, a minimum of 100 permutations are needed, so the threshold value is one that occurred one time in a hundred.

Below, I implement this procedure using the mcarlo function form the analogue package, using 1000 permutations.

Code
if(!"CordDmcarlo.rds"%in%dir("./Data/Novelty/")){
# use a Monte Carlo simulation of dissimilarities to define the suitability cut-off
  a<-Sys.time() # Clocks at ~4mins
  Maps0kaBP.CordDmcarlo <- mcarlo(Maps0kaBPTble, # current time Taxon data.frame 
                             method = "chord", # dissimilarity coefficient to use
                             nsamp = 1000, # number of permutations
                             type = "paired", # type of permutation or simulation to perform
                             replace = FALSE # sampling with replacement?
                             )
  saveRDS(Maps0kaBP.CordDmcarlo,"./Data/Novelty/CordDmcarlo.rds")
  a-Sys.time()
}
if(!"Maps0kaBP.CordDmcarlo"%in%ls()){
  Maps0kaBP.CordDmcarlo <- readRDS("./Data/Novelty/CordDmcarlo.rds")
}
Maps0kaBP.CordDmcarlo

    Simulated Dissimilarities

Simulation type : paired 
No. simulations : 1000 
Coefficient     : chord 

Summary of simulated distribution:
    Min 1st Qu.  Median    Mean 3rd Qu.     Max 
 0.0389  0.4358  1.0072  1.0189  1.5992  2.3579 

Percentiles of simulated distribution:
    1%   2.5%     5%    10%    90%    95%  97.5%    99% 
0.0492 0.0719 0.1196 0.1924 1.8323 1.9849 2.0811 2.1745 

Based on the estimation above, the cut-off value for non-analog growth form assemblages is 0.0492. Any dissimilarity above that value would indicate a non-analogue assemblage.

Now, using this value, the CordDistRast object, which contains the dissimilarity values could be masked to indicate which of the current areas are novel when compared to past conditions (Figure 5).

Code
# defining the suitable cut-off value
CutOffValCordD <- quantile(Maps0kaBP.CordDmcarlo,0.01)
# Plot the cut-off map
plot(CordDistRast > CutOffValCordD,
     main = "1950's Non-analogue areas compared to 21kaBP\n[Cord-Dist - Monte-Carlo based cut-off]")
plot(wrld_simpl, add=T)

Figure 5: Areas where NOVEL growth form compositional assemblages are expected based on a cord-distance/ROC-defined criteria.

If you would like to see the cord-distance/ROC-defined cutoff criteria in the context of all minimum cord-distances you can plot the histogram of distances (Figure 6):

Code
# Histogram of dissimilarities representing the distribution of dissimilarity distances and the cut-off values
hist(CordDistRast)
# Cut-off value
abline(v=CutOffValCordD)
legend("topright",
       paste0("cut-off = ",
              round(CutOffValCordD,4)))

Figure 6: Histogram showing the distibution of distances and the cord-distance/ROC-defined cutoff.

The same can now be done for the novelty estimated using the chi-squared distance. Below is the visualization for a ChiSrq-distance/MonteCarlo-defined cut-off criteria (Figure 7):

Code
if(!"CordDmcarlo.rds"%in%dir("./Data/Novelty/")){
# use a Monte Carlo simulation of dissimilarities to define the suitability cut-off
  a<-Sys.time() # Clocks at ~4mins
  Maps0kaBP.ChiSqmcarlo <- mcarlo(Maps0kaBPTble, # current time Taxon data.frame 
                             method = "chi.square", # dissimilarity coefficient to use
                             nsamp = 1000, # number of permutations
                             type = "paired", # type of permutation or simulation to perform
                             replace = FALSE # sampling with replacement?
                             )
  a-Sys.time()
  saveRDS(Maps0kaBP.ChiSqmcarlo,"./Data/Novelty/ChiSqmcarlo.rds")
}
if(!"Maps0kaBP.ChiSqmcarlo"%in%ls()){
  Maps0kaBP.ChiSqmcarlo <- readRDS("./Data/Novelty/ChiSqmcarlo.rds")
}

# defining the suitable cut-off value
CutOffValChiSqr <- quantile(Maps0kaBP.ChiSqmcarlo,0.01)
# Plot the cut-off map
plot(ChiDistRast > CutOffValChiSqr,
     main = "1950's Non-analogue areas compared to 21kaBP\n[Chi-Sqrd - Monte-Carlo based cut-off]")
plot(wrld_simpl, add=T)

Figure 7: Areas where NOVEL growth form compositional assemblages are expected based on a ChiSrq-distance/MonteCarlo-defined criteria.

Again, if you would like to see the ChiSrq-distance/MonteCarlo-defined cutoff criteria in the context of all minimum ChiSrq-distances you can plot the histogram of distances (Figure 8):

Code
## Histogram of dissimilarities representing the distribution of dissimilarity distances and the cut-off values
hist(ChiDistRast)
# cut-off value
abline(v=CutOffValChiSqr)
legend("topright",
       paste0("cut-off = ",
              round(CutOffValChiSqr,4)))

Figure 8: Histogram showing the distibution of distances and the ChiSrq-distance/MonteCarlo-defined cutoff.

An alternative form of defining the cut-off is to use the Receiver Operating Characteristic (ROC) curve, and I can divide the current conditions a priori into types of samples (e.g. vegetation types). Based don this, a site is an analogue for another site if they belong to the same group, and not an analogue if they come from different groups.

ROC curves are drawn using two measures of performance: i) sensitivity: the proportion of true analogues out of all sites said to be analogues on the basis of the cut-off - drawn on the y-axis. ii) specificity: the proportion of true non-analogues out of all non-analogues drawn on x-axis.

Here, like in species distribution modelling, the goal is to define a cut-off value that minimizes the false positive error (classifying two non-analogous samples as analogues) and the false negative error (classifying two analogous samples as non-analogues). That point is where mis-classifications are low: the True Positive Rate (i.e. sensitivity) are high, and Positive Rate (1-specificity) are low.

Below, I implement this procedure using the roc function form the analogue package, using 1000 permutations:

Code
if(!"CordDROC.rds"%in%dir("./Data/Novelty/")){
# load the classified map
  dbiome <- rast("./Data/biome_mclust_nodapc_18.tif")
# Nearest neighbor sample to the same extend as the growth form map
  dbiome <- resample(dbiome,
                     Maps0kaBP,
                     method = "near")
## Generate a vector of cells with class identity
  BiomeID <- values(dbiome,
                    dataframe = TRUE, # Make it a Data.frame
                    na.rm = FALSE # Keep NAs
                    )

# Estimate the ROC threshold
  a<-Sys.time() # Clocks ~5Min
  Map0k.CordDist <- analogue::distance(Maps0kaBPTble, method = "chord")
  Map0k.CordD.ROC <- roc(Map0k.CordDist, # current time Taxon data.frame 
                         groups = BiomeID[as.numeric(row.names(Maps0kaBPTble)),] # vector of group memberships
                         )
  a-Sys.time()
  saveRDS(Map0k.CordD.ROC,"./Data/Novelty/CordDROC.rds")
}
if(!"Map0k.CordD.ROC"%in%ls()){
  Map0k.CordD.ROC <- readRDS("./Data/Novelty/CordDROC.rds")
}

Based on the estimation above, the cut-off value for non-analog growth form assemblages is 0.0559. Any dissimilarity above that value would indicate a non-analogue assemblage.

Now, like with the Monte Carlo approach, I can classify the CordDistRast object as areas that are(are) novel when compared to past conditions (Figure 9).

Code
## defining the suitable cut-off value
CordDCutOffVal <- Map0k.CordD.ROC$statistics["Combined","Opt. Dis."]
# Plot the cut-off map
plot(CordDistRast > CordDCutOffVal,
     main = "1950's Non-analogue areas compared to 21kaBP\n[Cord Dist - ROC based cut-off]",
     xpd = NA)
plot(wrld_simpl, add=T)

Figure 9: Areas where NOVEL growth form compositional assemblages are expected based on a cord-distance/ROC-defined criteria.

I can also show the cutoff value in the context of all distances (Figure 10).

Code
# Histogram of dissimilarities representing the distribution of dissimilarity distances and the cut-off values
hist(CordDistRast)
# cut-off value
abline(v=CordDCutOffVal)
legend("topright",
       paste0("cut-off = ",
              round(CordDCutOffVal,4)))

Figure 10: Histogram showing the distibution of distances and the cord-distance/ROC-defined cutoff.

Now, the same procedure is done using a Chi-Squared distance (Figure 11).

Code
if(!"ChiDROC.rds"%in%dir("./Data/Novelty/")){
# Estimate the ROC threshold
  a<-Sys.time()  # Clocks ~5Min
  Map0k.ChiDist <- analogue::distance(Maps0kaBPTble, method = "chi.square")
  Map0k.ChiD.ROC <- roc(Map0k.ChiDist, # current time Taxon data.frame 
                        groups = BiomeID[as.numeric(row.names(Maps0kaBPTble)),] # vector of group memberships
                        )
  a-Sys.time()
  saveRDS(Map0k.ChiD.ROC,"./Data/Novelty/ChiDROC.rds")
}
if(!"Map0k.ChiD.ROC"%in%ls()){
  Map0k.ChiD.ROC <- readRDS("./Data/Novelty/ChiDROC.rds")
}
Map0k.ChiD.ROC

    ROC curve of dissimilarities

Discrimination for all groups:

Optimal Dissimilarity = 0.079 

AUC = NA, p-value: < 2.22e-16
No. within: 42663   No. outside: 725271 
Code
# defining the suitable cut-off value
ChiDCutOffVal <- Map0k.ChiD.ROC$statistics["Combined","Opt. Dis."]
# Plot the cut-off map
plot(ChiDistRast > ChiDCutOffVal,
     main = "1950's Non-analogue areas compared to 21kaBP\n[Chi Dist - ROC based cut-off]",
     xpd=NA)
plot(wrld_simpl, add=T)

Figure 11: Areas where NOVEL growth form compositional assemblages are expected based on a ChiSqr-distance/ROC-defined criteria.

Now, the same procedure is done using a cord-distance (Figure 12).

Code
# Histogram of dissimilarities representing the distribution of dissimilarity distances and the cut-off values
hist(ChiDistRast,
     main = "Compositional distance\n0kaBP to 21kaBP [Evergreen trees]",
     cex.main = 1.5)
# cut-off value
abline(v=ChiDCutOffVal)
legend("topright",
       paste0("cut-off = ",
              round(ChiDCutOffVal,4)))

Figure 12: Histogram showing the distibution of distances and the cord-distance/ROC-defined cutoff.

Note that when using the ROC defined threshold the maps done with Cord-distance and Chi-Squared distances are almost the same (there are differences but as a whole are negligible).

Velocity of phytoclimatic change as a metric of Ecological Novelty.

This mechanism focuses on measuring how fast would the “suitability” surface of for a given taxa would move in space - assuming that taxa would move from an area of low to an area of higher suitability between two times points. Estimations flow the implementation in (Ordonez et al. 2014) and (Ordonez, Williams, and Svenning 2016).

The approach used to estimate the Velocity of phytoclimatic change (that is the magnitude and direction of the change vector). Build on the approach developed by (Loarie et al. 2009), were velocity for an environmental variable (e.g., temperature) is estimated as:

\[V_{l} = \frac{\text{d}c/\text{d}t}{\text{d}c/\text{d}x}\]

where \(\frac{\text{d}c}{\text{d}t}\) is the the ration between the projected change per unit time; and \(\frac{\text{d}c}{\text{d}x}\) is the local spatial gradient in the variable of interest.

Here, I apply this approach to each growth form suitability maps rather than to a single climate variable (like in (Serra-Diaz et al. 2014)).

Estimating temporal gradients of change

I start by estimating the temporal gradient (i.e., \(\frac{\text{d}c}{\text{d}t}\) ) that represents projected change per unit time.

This can be estimated in different ways:

  1. As the slope of a generalized least squares regression on suitability maps for the 21kaBP to 0kaBP period for each cell with a autocorrelation structure of order one (AR1 model). This is the approach taken in (Dobrowski et al. 2013) and (Ordonez, Williams, and Svenning 2016).

  2. As the slope of a linear model on suitability maps for the 21kaBP to 0kaBP period for each cell. This is the approach taken in by (Loarie et al. 2009).

  3. The anomaly between the to end periods (21kaBP vs. 0kaBP) divided by the time between the start (21kaBP) and end (0kaBP) points. This is the approach take in (Sandel et al. 2011).

  4. A new approach we use here is based on the median differences between consecutive time periods. You can think of this as a rough way to consider the non linear trends in the relative changes in suitability over the evaluated time span.

  5. A second new approach we use here is based on the sum of absolute differences between consecutive time periods divided by the time between the start (21kaBP) and end (0kaBP) points. You can think of this as a rough way to consider the total amount of change over the evaluated time span.

A point to highlight is that the significance in the temporal trend can only be estimated for regression based approaches. For these, significance can be defined based on the p-value of the model regression slope but for now we will not take that in to consideration.

I developed the TempGradFnc() function to estimate these using all five models of temporal change. Outputs for these, can be seen HERE. For illustration, (Figure 13) shows how these can be estimated for all growth forms with the median differences between consecutive time periods approach.

Code
if(length(dir("./Data/Velocity/TempChng_Anomaly1/", pattern = "All_"))==0){
#Loop over all Growth Forms
  TempHetByGF <- lapply(NamesDtFrm$Acro2,
                        function(GF.Use){#(GF.Use <- NamesDtFrm$Acro2[1])
## **Zero**: Generate a SpatRaster with the suitability rasteres
        MapsList <- lapply(1:44,
                           function(i){
                             rasterize(x = as.matrix(xy), # points as a matrix
                                       y = cr.ea, #  template raster
                                       values = pml[[i]][,GF.Use] # Values to aggregate
                             )
                           })
# make a multi-band SpatRaster
        MapsPer.GF <- do.call("c",MapsList)
## **First**: Generate a SpatRaster that estimates the temporal trend for each cell using the app function from terra
        TempHetRast <- app(MapsPer.GF,
                           fun = function(i, ff) ff(i,method="Anomaly1"),#
                           ff=TempGradFnc,
                           cores = 10 # the function is run automatically in parallel in 10 cores
        )
        TempHetRast <- mask(TempHetRast, # Mask based on the input raster
                            MapsPer.GF[[1]])
        return(TempHetRast)
        })
} else{
  TempHetByGF <- lapply(NamesDtFrm$Acro2,
                        function(GF.Use){#(GF.Use <- NamesDtFrm$Acro2[1])
                          rast(paste0("./Data/Velocity/TempChng_lm/All_",
                                      GF.Use,
                                      "_TempChng.tif"))  
                        })
  }
# Merge all values 
TempHetByGF <- do.call("c",TempHetByGF)
names(TempHetByGF) <- NamesDtFrm$Name

# plot the Temporal gradient
levelplot(raster::stack(TempHetByGF), # level plot only works with raster files
          scales = list(draw=FALSE), # To remove the Latlong
          main=list("Temporal gradient [% per 100 yrs]\n(median difference of inter period diferecnes)",side=1,line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(col=rev(hcl.colors(100, palette = "RdYlBu"))) # add a legend
          ) +
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

Figure 13: Temporal heterogeneity (%-change per 100yrs) between 21kaBA and 0kaBP

To get a sense of how these changes compare between methods of evaluating temporal heterogeneity, I now show all five estimates for Evergreen trees (Figure 14).

Code
# Load the Temp Heterogeneity surfaces
TempHetTE <- lapply(dir("./Data/Velocity",pattern = "TempChng"),
                        function(TimePer){#(TimePer <- dir("./Data/Velocity",pattern = "TempChng")[1])
                          rast(paste0("./Data/Velocity/",
                                      TimePer,
                                      "/All_TE_TempChng.tif"))  
                        })
# Merge all values 
TempHetTE <- do.call("c",TempHetTE)
names(TempHetTE) <- gsub("TempChng_","",dir("./Data/Velocity",pattern = "TempChng"))

# plot the Temporal gradient
levelplot(raster::stack(abs(TempHetTE)^(1/5)), # level plot only works with raster files
          scales = list(draw=FALSE), # To remove the Latlong
          main=list("Absolute Temporal gradient [% per 100 yrs]\n(Estimated using five aproaches - Sacled as $X^(1/5)$ for ease of visualization)",side=1,line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(col=rev(hcl.colors(100, palette = "RdYlBu"))) # add a legend
          ) +
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

Figure 14: Contrast of ABSOLUTE temporal heterogeneity (%-change per 100yrs) metrics for Evergreen trees. Scaled as \(X^ rac{1}{5}\) for ease of visualization.

Estimating the spatial heterogeneity

The second step is to estimating the spatial heterogeneity (change per unit space; i.e., \(\frac{\text{d}c}{\text{d}x}\)) as in (Burrows et al. 2011), (Dobrowski et al. 2013) and (Ordonez, Williams, and Svenning 2016) for each map cell as “the slope of proportions” using the maximum average technique ((Burrough, McDonnell, and Lloyd 2015)).

The method focuses on estimating the average change in the West-East (W-E) direction (negative values indicate a westward direction), and the North-South (N-S) direction (negative values indicate a equatorial direction) and divided by the avg distance between the cells (47km in the West-East direction and 66km in the North-South direction). The overall spatial gradients is then calculated as the vector sum of the N-S and W-E gradients, with the associated vector angle giving the direction of the gradient (Figure 15).

Code
if(length(dir("./Data/Velocity/SpatHet/",
              pattern="All_"))==0){
#Loop over all Growth Forms
  SpaceHetByGF <- lapply(NamesDtFrm$Acro2,
                        function(GF.Use){#(GF.Use <- NamesDtFrm$Acro2[1])
## **Zero**: Generate a SpatRaster with the suitability rasteres
        MapsList <- lapply(1:44,
                           function(i){
                             rasterize(x = as.matrix(xy), # points as a matrix
                                       y = cr.ea, #  template raster
                                       values = pml[[i]][,GF.Use] # Values to aggregate
                             )
                           })
# make a multi-band SpatRaster
        MapsPer.GF <- do.call("c",MapsList)
## **Second**: Estimate the spatial gradients magnitude using a using the maximum average technique [Burrough & McDonnell 1998].
      SpaceHetRast <- SpatHetFnc(MapsPer.GF[[1]])
       return(SpaceHetRast)
        })
  } else{
    SpaceHetByGF <- lapply(NamesDtFrm$Acro2,
                        function(GF.Use){#(GF.Use <- NamesDtFrm$Acro2[1])
                          rast(paste0("./Data/Velocity/SpatHet/All_",
                                      GF.Use,
                                      "_SpatHet.tif"))  
                        })
  }      
# Merge all values 
SpaceHetByGF <- do.call("c",SpaceHetByGF)
names(SpaceHetByGF) <- NamesDtFrm$Name

# plot the Spatial gradient
levelplot(raster::stack((SpaceHetByGF)^(1/5)), # level plot only works with raster files
          scales = list(draw=FALSE), # To remove the Latlong
          main=list("Spatial gradient [(% per kg)^(1/5)]",side=1,line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(col=rev(hcl.colors(100, palette = "RdYlBu"))) # add a legend
          ) +
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

Figure 15: Spatial heterogeneity (%-change per km) for the baseline period (21kaBP). Values Scaled as \(X^ rac{1}{5}\) to increse contrast.

Estimating the magnitude of the velocity vectors.

With these two variables (that is the spatial heterogeneity, and the temporal heterogeneity) I can estimate the velocity as the ratio between these two. Here there is a need to do some corrections based on some mathematical issues of estimating ratios.

One of this is estimating a ration when the denominator is zero (a case here when the spatial gradient is zero as all neighboring cells have the same value). In this situation you can either define these cells to have the minimum value of differentiation allowed based on the “accuracy” of your measurements (as done by by (Loarie et al. 2009) and (Sandel et al. 2011)). An alternative, is to define these areas as the “maximum” velocity used for visualization (1000km per 100yr in this case).

A second issue is velocities that are two small as the spatial change is considerably larger than the temporal change. Based on the spatiotemporal scale of the study, these changes can be set to the minimum possible velocity (0.001km per 100yr in this case).

With these clarifications, we then estimated the velocity of change for each of the evaluated growth forms (Figure 16). But as I also implemented multiple metrics of temporal change, velocity estimates need to be estimated for each of these.

Code
if(length(dir("./Data/Velocity/Velocity_Anomaly1/",
              pattern="All_"))==0){
#Loop over all Growth Forms
  VelocityByGF <- lapply(NamesDtFrm$Acro2,
                        function(GF.Use){#(GF.Use <- NamesDtFrm$Acro2[1])
## **Zero**: load a SpatRaster with the Temporal heterogeneity 
                          TempHetRast <-  rast(paste0("./Data/Velocity/Velocity_Anomaly1/All_",GF.Use,"_TempChng.tif"))
## **Zero**: load a SpatRaster with the Space heterogeneity 
                          SpaceHetRast <- rast(paste0("./Data/Velocity/SpatHet/All_",GF.Use,"_SpatHet.tif"))
## **Forth**: Estimate the velocity magnitude (i.e., speed) as the ratio between the temporal and spatial gradient.
                          Velocity <- VelocityFnc (TempHetRast,
                                                   SpaceHetRast)
                          return(Velocity)
                          })
  } else{
    VelocityByGF <- lapply(NamesDtFrm$Acro2,
                        function(GF.Use){#(GF.Use <- NamesDtFrm$Acro2[1])
                          rast(paste0("./Data/Velocity/Velocity_Anomaly1/All_",
                                      GF.Use,
                                      "_Velocity.tif"))
                        })
  } 
# Merge all values 
VelocityByGF <- do.call("c",VelocityByGF)
names(VelocityByGF) <- NamesDtFrm$Name

# plot the Temporal gradient
levelplot(raster::stack(log10(VelocityByGF)), # level plot only works with raster files
          at=seq(-3,3,length.out=100), # zlim 
          scales = list(draw=FALSE), # To remove the Latlong
          main=list("Velocity of climate changt [km / 100 yrs]",side=1,line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(labels = list(at = -3:3, # add a legend and its properties
                                        labels = 10^c(-3:3),
                                        col=rev(hcl.colors(100, palette = "RdYlBu"))))
          ) +
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

Figure 16: Velocity of change for Per Growth Form

To get a sense of how these changes compare between methods of evaluating temporal heterogeneity, I now show all five estimates for Evergreen trees (Figure 17).

Code
# Load the Temp Heterogeneity surfaces
VelocityTE <- lapply(dir("./Data/Velocity",pattern = "Velocity"),
                        function(TimePer){#(TimePer <- dir("./Data/Velocity",pattern = "TempChng")[1])
                          rast(paste0("./Data/Velocity/",
                                      TimePer,
                                      "/All_TE_Velocity.tif"))  
                        })
# Merge all values 
VelocityTE <- do.call("c",VelocityTE)
names(VelocityTE) <- gsub("Velocity_","",dir("./Data/Velocity",pattern = "Velocity"))

# plot the Temporal gradient
levelplot(raster::stack(log10(VelocityTE)), # level plot only works with raster files
          at=seq(-3,3,length.out=100), # zlim
          scales = list(draw=FALSE), # To remove the Latlong
          main=list("Absolute Temporal gradient [% per 100 yrs]\n(Estimated using five aproaches)",side=1,line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(labels = list(at = -3:3,# add a legend and its properties
                                      labels = 10^c(-3:3),
                                      col=rev(hcl.colors(100, palette = "RdYlBu"))))
          ) +
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

Figure 17: Velocity of change (km per 100yrs) between 21kaBP and 0kaBP for Evergreen trees.

Estimating the Bearing of the velocity vector

Velocity as a vector has two components, a magnitude (defined here by the ration between temporal and spatial changes) and a direction (defined by the spatial gradient vector and the direction of the temporal change).

The calculation of the velocity vector bearing is determined by the East-West and North-South components of the spatial gradient. The sum of these define the magnitude and direction of vector of the spatial change (the one defining towards where the incline would move). The bearing of this resulting vector is the angle measured clockwise with 90o centered on the corresponding pole. Therefore, an angle ranging between 0 and 180 means the vector is Polewards bound and angles between 180 and 359 degrees is Equatorward bound. A 0 bearing means the vector is Eastbound in both the Northern and Southern hemisphere, and an angle of 180 is Westbound both the Northern and Southern hemisphere.

Furthermore, the direction of change is considered to be from high-values to low-values. So if the temporal change indicates that a cell is decreasing in values over time, the bearing of the vector would need to be filliped as it will change from being a source to being a sink.

I estimated the bearing of the velocity vector using the vector resulting from the sum of the rise and run run of the spatial gradient, reversing the direction of change is the temporal trend was negative (Figure 18).

Code
if(length(dir("./Data/Velocity/SpatHet/",
              pattern="All_"))==0){
#Loop over all Growth Forms
  SpaceHetByGF <- lapply(NamesDtFrm$Acro2,
                        function(GF.Use){#(GF.Use <- NamesDtFrm$Acro2[1])
## **Zero**: Generate a SpatRaster with the suitability rasteres
        MapsList <- lapply(1:44,
                           function(i){
                             rasterize(x = as.matrix(xy), # points as a matrix
                                       y = cr.ea, #  template raster
                                       values = pml[[i]][,GF.Use] # Values to aggregate
                             )
                           })
# make a multi-band SpatRaster
        MapsPer.GF <- do.call("c",MapsList)
## **Third**: Estimate the spatial gradients direction using a using the maximum average technique [Burrough & McDonnell 1998].
##            Here 90 degrees is poleward direction, so that 0 degrees is East in the north and West in the south
      a<-Sys.time()
      BearingByGF <- BearingFnc(MapsPer.GF)
       return(BearingByGF)
        })
  } else{
    BearingByGF <- lapply(NamesDtFrm$Acro2,
                        function(GF.Use){#(GF.Use <- NamesDtFrm$Acro2[1])
                          rast(paste0("./Data/Velocity/Bearing/All_",
                                      GF.Use,
                                      "_Bearing.tif"))  
                        })
  }      
# Merge all values 
BearingByGF <- do.call("c",BearingByGF)
names(BearingByGF) <- NamesDtFrm$Name      
      
# plot the Spatial gradient
levelplot(raster::stack(BearingByGF), # level plot only works with raster files
          scales = list(draw=FALSE), # To remove the Latlong
          main=list("Bearing [Degrees from North]",side=1,line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(at=seq(0,360,length.out=100),  # add a legend and its properties
                        labels = list(at = seq(0,360,by=90),
                                      labels = c("Eastward", "Poleward","Westward","Eastward"),
                                      col=rev(hcl.colors(100, palette = "RdYlBu"))))
          ) +
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

Figure 18: Bearing of the spatial gradients

Displacement of phytoclimatic change as a metric of Ecological Novelty.

As defined in (Ordonez, Williams, and Svenning 2016), displacement of climatic vectors is the mean of multiple velocity vectors. This metric indicates how “fast” would/will an environmental setup (weighting all variables equally) would move given the balance of local environmental gradients and temporal trends in environmental variables. A fast displacement suggests that the magnitudes of velocity vectors (i.e., speed) are rather large, whereas slow displacement indicate that the magnitudes of velocity vectors is small across evaluated growth forms.

Here, I leverage the velocity estimates for the 14 evaluated growth forms to assess how fast has the “growth form setup” changed since the Last Glacial Maxim (LGM; ~21ka) (?@fig-DisplacementAll). The Displacement of growth form suitability. For this, I used the geometric mean of the velocity vectors (here implemented as the mean of Log-10 velocities).

Code
if(length(dir("./Data/Displacement/Displacement_AbsDif",
              pattern="All_"))==0){
    Displacement <- lapply(dir("./Data/Displacement"),
                          function(Model){#(Model <- dir("./Data/Displacement")[1])
                            VelocityByGF <- lapply(NamesDtFrm$Acro2,
                                                   function(GF.Use){#(GF.Use <- NamesDtFrm$Acro2[1])
                                                     rast(paste0("./Data/Velocity/",
                                                                 gsub("Displacement_",
                                                                      "Velocity_",
                                                                      Model),
                                                                 "/All_",
                                                                 GF.Use,
                                                                 "_Velocity.tif"))
                                                   })
# Merge all values 
                            VelocityByGF <- do.call("c",VelocityByGF)
                            names(VelocityByGF) <- NamesDtFrm$Name
# Estimate the Displacement <-  geometric mean of velocities
                            Displacement <- 10^mean(log10(VelocityByGF))
                            return(Displacement)
                            })
    } else{
# Load the Displacement vectors
  Displacement <- lapply(dir("./Data/Displacement"),
                          function(Model){#(Model <- dir("./Data/Displacement")[1])
                            rast(paste0("./Data/Displacement/",Model,"/All_Displacement.tif"))
                          })
  }
# Merge all values 
Displacement <- do.call("c",Displacement)

# Give names to the layers
names(Displacement) <- gsub("Displacement_",
                            "",
                            dir("./Data/Displacement"))  

# plot the Displacement gradient
levelplot(raster::stack(log10(Displacement)), # level plot only works with raster files
          scales = list(draw=FALSE), # To remove the Latlong
          main=list("Displacement as the average over GF Velocities [km/ yrs]\n(Estimated using five aproaches)",side=1,line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(at=seq(-3,3,length.out=100),  # add a legend and its properties
                        labels = list(at = -3:3,
                                      labels = 10^c(-3:3),
                                      col=rev(hcl.colors(100, palette = "RdYlBu"))))
          ) +
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

Displacement of growth form suitability velocity vectors

Divergence of phytoclimatic velocity vecrtors as a metric of Ecological Novelty.

As defined in (Ordonez, Williams, and Svenning 2016), divergence of climatic vectors is the angle between two velocity vectors. In a multivariate context this definition is not practical, and I instead use the standard deviations of bearings as in (Burke et al. 2019). This metric indicates how “variable” are the directions of velocity vectors in a given location for the 14 evaluated growth forms. A low divergence suggests that all velocity vectors of all climate variables would tend to shift along the same axis of direction, whereas high divergences indicate that velocity vectors lack congruence in orientation and hence so might species distribution shifts.

Here, I leverage the velocity estimates for the 14 evaluated growth forms to assess how fast has the “growth form setup” changed since the Last Glacial Maxim (LGM; ~21ka) (Figure 19). The Displacement of growth form suitability. For this, I used the geometric mean of the velocity vectors (here implemented as the mean of Log-10 velocities).

Code
if(!"All_Divergence.tif"%in%dir("./Data/Divergence")){
  BearingByGF <- lapply(NamesDtFrm$Acro2,
                           function(GF.Use){#(GF.Use <- NamesDtFrm$Acro2[1])
                             rast(paste0("./Data/Velocity/Bearing/All_",
                                         GF.Use,
                                         "_Bearing.tif"))
                           })
# Estimate the Displacement <-  geometric mean of velocities
    Divergence <- app(BearingByGF,sd,na.rm=T)
} else{
  Divergence <- rast("./Data/Divergence/All_Divergence.tif")
}   
#plot the Divergence
plot(Divergence,
     main = "Divergence accorss growth form suitability",
     col = rev(hcl.colors(100,"RdYlBu")),
     xpd=NA)
# add world map
plot(wrld_simpl, add=T)
# add the ice maps 
plot(Ice, 
    col = gray(0.3,alpha = 0.6),
    legend = F,
      add = T)

Figure 19: Divergence of growth form suitability velocity vectors

Velocity per Greenland Interstadial and Greenland Stadials

Here I will estimate the velocity of change of the growth form suitability surfaces for specific time periods for the Last Termination defined by the GRIP Greenland ice-core oxygen isotope signal as in (Björck et al. 1998).

The four (4) periods to be used here correspond to two stadial episodes (Greenland Stadials 1 (GS-1) and 2 (GS-2)), one interstadial event (Greenland Interstadial 1 (GI-1)), and the Holocene. The GI-1 and GS-2 correspond to other known periods, namely the Bølling–Allerød warming (~ 14.7kaBP to 12.7kaBP) that corresponds to GI-1, and the subsequent Younger Dryas cooling (~ 12.7kaBP to 11.5kaBP) that corresponds to GS-1.

Table 1 - Depths and preliminary ages (ad1950) of the onset of events and episodes in the GRIP ice-core, including the Holocene epoch. Dates based on (Björck et al. 1998).

Events Ice-core age
Holocene 11.5kaBP to 0 kaBP
GS-1 12.7kaBP to 11.5kaBP
GI-1 14.7kaBP to 12.7kaBP
GS-2 21.2kaBP to 14.7kaBP

Like above, I also assess the velocity of change as the ratio between the spatial and temporal changes, for each growth form each time period. For this, I loop over the the time periods and growth forms, and divide the process in three stages, and in each stage. The codes is as the one in the Estimating the magnitude of the velocity vectors. Below I show the different velocity metrics (based on the median of between periods changes) across growth forms for each period (Figure 20; Figure 21, Figure 20, Figure 23):

Code
# label: VelocityPerPeriod
# Table with all the dates
# I will loop over all periods 
VelocityAllList <- lapply( c("Holocene","GS1","GI1","GS2"),
                          function(TimePer){#(TimePer<-"Holocene")
# Load and plot the Velocity magnitudes as a list
  VelocityAll <- lapply(NamesDtFrm$Acro2,
                        function(GF.Use){
                          rast(paste0("./Data/Velocity/Velocity_Anomaly1/",
                                      TimePer,"_",
                                      GF.Use,
                                      "_Velocity.tif"))
                          })
# Make a multi-band SpatRaster
  VelocityAll <- do.call("c",VelocityAll)
  names(VelocityAll) <- NamesDtFrm$Name
  return(VelocityAll)
})
names(VelocityAllList) <- c("Holocene","GS1","GI1","GS2")
Code
# plot the Velocity
levelplot(raster::stack(log10(VelocityAllList[["GS2"]])), # level plot only works with raster files
          scales = list(draw=FALSE), # To remove the Latlong
          at=seq(-3,3.000434,length.out=100), 
          main=list("Velocity [km per 100yrs]\nGS2",
                    side=1,
                    line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(labels = list(at = -3:3, # add a legend and its properties
                                      labels = 10^c(-3:3),
                                      col=rev(hcl.colors(100, palette = "RdYlBu"))))
          ) + 
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

Figure 20: Velocity for all evaluated growth forms during the Greenland Stadial-2 (Post-LGM)

Code
# plot the Velocity
levelplot(raster::stack(log10(VelocityAllList[["GI1"]])), # level plot only works with raster files
          at = seq(-3,3.1,length.out=100),
          scales = list(draw=FALSE), # To remove the Latlong
          main=list("Velocity [km per 100yrs]\nGI1",
                    side=1,
                    line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(labels = list(at = -3:3, # add a legend and its properties
                                      labels = 10^c(-3:3),
                                      col=rev(hcl.colors(100, palette = "RdYlBu"))))
          ) + 
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

Figure 21: Velocity for all evaluated growth forms during the Greenland Inerstadial-1 (BA)

Code
# plot the Velocity
levelplot(raster::stack(log10(VelocityAllList[["GS1"]])), # level plot only works with raster files
          scales = list(draw=FALSE), # To remove the Latlong
          at=seq(-3,3.000434,length.out=100), 
          main=list("Velocity [km per 100yrs]\nGS1",
                    side=1,
                    line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(labels = list(at = -3:3, # add a legend and its properties
                                      labels = 10^c(-3:3),
                                      col=rev(hcl.colors(100, palette = "RdYlBu"))))
          ) + 
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

Figure 22: Velocity for all evaluated growth forms during the Greenland Stadial-1 (YD)

Code
# plot the Velocity
levelplot(raster::stack(log10(VelocityAllList[["Holocene"]])), # level plot only works with raster files
          scales = list(draw=FALSE), # To remove the Latlong
          at=seq(-3,3.000434,length.out=100), 
          main=list("Velocity [km per 100yrs]\nHolocene",
                    side=1,
                    line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(labels = list(at = -3:3, # add a legend and its properties
                                      labels = 10^c(-3:3),
                                      col=rev(hcl.colors(100, palette = "RdYlBu"))))
          ) + 
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

Figure 23: Velocity for all evaluated growth forms during the Holocene (YD)

Displacement and divergence of phytoclimatic change per Greenland Interstadial and Greenland Stadials.

Now with the metrics of velocity and bearings I can then estimate displacement and divergence of climatic vectors (cf. (Ordonez, Williams, and Svenning 2016)) for mutiple periods (Figure 24, Figure 25, Figure 26, (Figure 27).

As a reminder, displacement of climatic vectors is the mean of multiple velocity vectors and indicates how “fast” would/will an environmental setup (weighting all variables equally) would move given the balance of local environmental gradients and temporal trends in environmental variables.

Code
(TimePer<-"GS2")
[1] "GS2"
Code
# Load and plot the Velocity magnitudes as a list
  DisplacementAllList <- lapply(c("AbsDif", "Anomaly1", "Anomaly2", "ARM", "lm"),
                           function(Model){#(Model<-"AbsDif")
                             rast(paste0("./Data/Displacement/Displacement_",
                                         Model,
                                         "/",
                                         TimePer,
                                         "_Displacement.tif"))
                           })
# Make a multi-band SpatRaster
  DisplacementAll <- do.call("c",DisplacementAllList)
  names(DisplacementAll) <- c("AbsDif", "Anomaly1", "Anomaly2", "ARM", "lm")
# plot the Displacement
  print(
  levelplot(raster::stack(log10(DisplacementAll)), # level plot only works with raster files
            at=seq(-3,3.000434,length.out=100),
            scales = list(draw=FALSE), # To remove the Latlong
            main=list(paste0("Displacement [km per 100yrs]\n",TimePer),
                      side=1,
                      line=-0.5), # Main title
            col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
            colorkey = list(labels = list(at = -3:3, # add a legend and its properties
                                          labels = 10^c(-3:3),
                                          col=rev(hcl.colors(100, palette = "RdYlBu"))))
            ) + 
# Add the country outlines
  layer(sp.lines(wrld_simpl)
        ) + 
# Plot the Ice maps 
  levelplot(raster::raster(Ice), # add the ice maps 
            col.regions = "dark grey",
            alpha.regions = 0.6)
)

Figure 24: *Displacement of growth form suitability velocity vectors during the GS2**

Code
(TimePer<-"GI1")
[1] "GI1"
Code
# Load and plot the Velocity magnitudes as a list
  DisplacementAllList <- lapply(c("AbsDif", "Anomaly1", "Anomaly2", "ARM", "lm"),
                           function(Model){#(Model<-"AbsDif")
                             rast(paste0("./Data/Displacement/Displacement_",
                                         Model,
                                         "/",
                                         TimePer,
                                         "_Displacement.tif"))
                           })
# Make a multi-band SpatRaster
  DisplacementAll <- do.call("c",DisplacementAllList)
  names(DisplacementAll) <- c("AbsDif", "Anomaly1", "Anomaly2", "ARM", "lm")
# plot the Displacement
  print(
  levelplot(raster::stack(log10(DisplacementAll)), # level plot only works with raster files
            at=seq(-3,3.000434,length.out=100),
            scales = list(draw=FALSE), # To remove the Latlong
            main=list(paste0("Displacement [km per 100yrs]\n",TimePer),
                      side=1,
                      line=-0.5), # Main title
            col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
            colorkey = list(labels = list(at = -3:3, # add a legend and its properties
                                          labels = 10^c(-3:3),
                                          col=rev(hcl.colors(100, palette = "RdYlBu"))))
            ) + 
# Add the country outlines
  layer(sp.lines(wrld_simpl)
        ) + 
# Plot the Ice maps 
  levelplot(raster::raster(Ice), # add the ice maps 
            col.regions = "dark grey",
            alpha.regions = 0.6)
)

Figure 25: *Displacement of growth form suitability velocity vectors during the GI1**

Code
(TimePer<-"GS1")
[1] "GS1"
Code
# Load and plot the Velocity magnitudes as a list
  DisplacementAllList <- lapply(c("AbsDif", "Anomaly1", "Anomaly2", "ARM", "lm"),
                           function(Model){#(Model<-"AbsDif")
                             rast(paste0("./Data/Displacement/Displacement_",
                                         Model,
                                         "/",
                                         TimePer,
                                         "_Displacement.tif"))
                           })
# Make a multi-band SpatRaster
  DisplacementAll <- do.call("c",DisplacementAllList)
  names(DisplacementAll) <- c("AbsDif", "Anomaly1", "Anomaly2", "ARM", "lm")
# plot the Displacement
  print(
levelplot(raster::stack(log10(DisplacementAll)), # level plot only works with raster files
            at=seq(-3,3.000434,length.out=100),
            scales = list(draw=FALSE), # To remove the Latlong
            main=list(paste0("Displacement [km per 100yrs]\n",TimePer),
                      side=1,
                      line=-0.5), # Main title
            col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
            colorkey = list(labels = list(at = -3:3, # add a legend and its properties
                                          labels = 10^c(-3:3),
                                          col=rev(hcl.colors(100, palette = "RdYlBu"))))
            ) + 
# Add the country outlines
  layer(sp.lines(wrld_simpl)
        ) + 
# Plot the Ice maps 
  levelplot(raster::raster(Ice), # add the ice maps 
            col.regions = "dark grey",
            alpha.regions = 0.6)
)

Figure 26: *Displacement of growth form suitability velocity vectors during the GS1**

Code
(TimePer<-"Holocene")
[1] "Holocene"
Code
# Load and plot the Velocity magnitudes as a list
  DisplacementAllList <- lapply(c("AbsDif", "Anomaly1", "Anomaly2", "ARM", "lm"),
                           function(Model){#(Model<-"AbsDif")
                             rast(paste0("./Data/Displacement/Displacement_",
                                         Model,
                                         "/",
                                         TimePer,
                                         "_Displacement.tif"))
                           })
# Make a multi-band SpatRaster
  DisplacementAll <- do.call("c",DisplacementAllList)
  names(DisplacementAll) <- c("AbsDif", "Anomaly1", "Anomaly2", "ARM", "lm")
# plot the Displacement
  print(
levelplot(raster::stack(log10(DisplacementAll)), # level plot only works with raster files
            at=seq(-3,3.000434,length.out=100),
            scales = list(draw=FALSE), # To remove the Latlong
            main=list(paste0("Displacement [km per 100yrs]\n",TimePer),
                      side=1,
                      line=-0.5), # Main title
            col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
            colorkey = list(labels = list(at = -3:3, # add a legend and its properties
                                          labels = 10^c(-3:3),
                                          col=rev(hcl.colors(100, palette = "RdYlBu"))))
            ) + 
# Add the country outlines
  layer(sp.lines(wrld_simpl)
        ) + 
# Plot the Ice maps 
  levelplot(raster::raster(Ice), # add the ice maps 
            col.regions = "dark grey",
            alpha.regions = 0.6)
)

Figure 27: *Displacement of growth form suitability velocity vectors during the Holocene**

Divergence of climatic vectors (?@fig-DivergencePerPeriod) is the standard deviations of bearings as in (Burke et al. 2019). This metric indicates how “variable” are the directions of velocity vectors in a given location for the 14 evaluated growth forms.

Code
# I will loop over all periods 
# Load and plot the Divergence magnitudes as a list
Divergence <- lapply(c("Holocene","GS1","GI1","GS2" ),
                         function(TimePer){
                           rast(paste0("./Data/Divergence/",
                                       TimePer,
                                       "_Divergence.tif"))
                         })
# Make a multi-band SpatRaster
Divergence <- do.call("c",Divergence)
names(Divergence) <- c("Holocene","GS1","GI1","GS2" )

# plot the Displacement
levelplot(raster::stack(Divergence), # level plot only works with raster files
          scales = list(draw=FALSE), # To remove the Latlong
          at = seq(0,180,length.out=100),
          main=list("Divergence of growth form suitability",
                    side=1,
                    line=-0.5), # Main title
          col.regions = rev(hcl.colors(100,"RdYlBu")), # set the colors
          colorkey = list(at=seq(0,180,length.out=100),  # add a legend and its properties
                        labels = list(at = seq(0,180,by=30),
                                      labels = seq(0,180,by=30),
                                      col=rev(hcl.colors(100, palette = "RdYlBu"))))
          ) + 
# Add the country outlines
layer(sp.lines(wrld_simpl)
      ) + 
# Plot the Ice maps 
levelplot(raster::raster(Ice), # add the ice maps 
          col.regions = "dark grey",
          alpha.regions = 0.6)

Divergence of growth form suitability velocity vectors during four Post-LGM periods

References

Björck, Svante, Michael J. C. Walker, Les C. Cwynar, Sigfus Johnsen, Karen-Luise Knudsen, J. John Lowe, Barbara Wohlfarth, and INTIMATE Members. 1998. “An Event Stratigraphy for the Last Termination in the North Atlantic Region Based on the Greenland Ice-Core Record: A Proposal by the INTIMATE Group.” Journal of Quaternary Science 13 (4): 283–92. https://doi.org/10.1002/(SICI)1099-1417(199807/08)13:4<283::AID-JQS386>3.0.CO;2-A.
Burke, Kevin D., John W. Williams, Simon Brewer, Walter Finsinger, Thomas Giesecke, David J. Lorenz, and Alejandro Ordonez. 2019. “Differing Climatic Mechanisms Control Transient and Accumulated Vegetation Novelty in Europe and Eastern North America.” Philosophical Transactions of the Royal Society B: Biological Sciences 374 (1788): 20190218. https://doi.org/10.1098/rstb.2019.0218.
Burrough, P. A., Rachael McDonnell, and Christopher D. Lloyd. 2015. Principles of Geographical Information Systems. Third edition. Oxford ; New York: Oxford University Press.
Burrows, Michael T., David S. Schoeman, Lauren B. Buckley, Pippa Moore, Elvira S. Poloczanska, Keith M. Brander, Chris Brown, et al. 2011. “The Pace of Shifting Climate in Marine and Terrestrial Ecosystems.” Science 334 (6056): 652–55. https://doi.org/10.1126/science.1210288.
Conradi, Timo, Jasper A. Slingsby, Guy F. Midgley, Henning Nottebrock, Andreas H. Schweiger, and Steven I. Higgins. 2020. “An Operational Definition of the Biome for Global Change Research.” New Phytologist 227 (5): 1294–1306. https://doi.org/10.1111/nph.16580.
Dobrowski, Solomon Z., John Abatzoglou, Alan K. Swanson, Jonathan A. Greenberg, Alison R. Mynsberge, Zachary A. Holden, and Michael K. Schwartz. 2013. “The Climate Velocity of the Contiguous United States During the 20th Century.” Global Change Biology 19 (1): 241–51. https://doi.org/10.1111/gcb.12026.
Flower, R. J., S. Juggins, and R. W. Battarbee. 1997. “Matching Diatom Assemblages in Lake Sediment Cores and Modern Surface Sediment Samples: The Implications for Lake Conservation and Restoration with Special Reference to Acidified Systems.” Hydrobiologia 344: 27–40. https://doi.org/10.1023/A:1002941908602.
Gavin, Daniel G, W.Wyatt Oswald, Eugene R Wahl, and John W Williams. 2003. “A Statistical Approach to Evaluating Distance Metrics and Analog Assignments for Pollen Records.” Quaternary Research 60 (3): 356–67. https://doi.org/10.1016/S0033-5894(03)00088-7.
Loarie, Scott R., Philip B. Duffy, Healy Hamilton, Gregory P. Asner, Christopher B. Field, and David D. Ackerly. 2009. “The Velocity of Climate Change.” Nature 462 (7276): 1052–55. https://doi.org/10.1038/nature08649.
Ordonez, Alejandro, Sebastián Martinuzzi, Volker C. Radeloff, and John W. Williams. 2014. “Combined Speeds of Climate and Land-Use Change of the Conterminous US Until 2050.” Nature Climate Change 4 (9): 811–16. https://doi.org/10.1038/nclimate2337.
Ordonez, Alejandro, John W. Williams, and Jens-Christian Svenning. 2016. “Mapping Climatic Mechanisms Likely to Favour the Emergence of Novel Communities.” Nature Climate Change 6 (12): 1104–9. https://doi.org/10.1038/nclimate3127.
Overpeck, J. T., T. Webb, and I. C. Prentice. 1985. “Quantitative Interpretation of Fossil Pollen Spectra: Dissimilarity Coefficients and the Method of Modern Analogs.” Quaternary Research 23 (1): 87–108. https://doi.org/10.1016/0033-5894(85)90074-2.
Sandel, B., L. Arge, B. Dalsgaard, R. G. Davies, K. J. Gaston, W. J. Sutherland, and J.-C. Svenning. 2011. “The Influence of Late Quaternary Climate-Change Velocity on Species Endemism.” Science 334 (6056): 660–64. https://doi.org/10.1126/science.1210173.
Serra-Diaz, Josep M., Janet Franklin, Miquel Ninyerola, Frank W. Davis, Alexandra D. Syphard, Helen M. Regan, and Makihiko Ikegami. 2014. “Bioclimatic Velocity: The Pace of Species Exposure to Climate Change.” Edited by Matt Fitzpatrick. Diversity and Distributions 20 (2): 169–80. https://doi.org/10.1111/ddi.12131.
Simpson, Gavin L. 2007. “Analogue Methods in Palaeoecology: Using the Analogue Package.” Journal of Statistical Software 22 (2). https://doi.org/10.18637/jss.v022.i02.
Williams, John W., and Stephen T. Jackson. 2007. “Novel Climates, No-Analog Communities, and Ecological Surprises.” Frontiers in Ecology and the Environment 5 (9): 475–82. https://doi.org/10.1890/070037.
Williams, John W., Stephen T. Jackson, and John E. Kutzbach. 2007. “Projected Distributions of Novel and Disappearing Climates by 2100 AD.” Proceedings of the National Academy of Sciences 104 (14): 5738–42. https://doi.org/10.1073/pnas.0606292104.